setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
library(dplyr)
library(ggplot2)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
id1=read.table("../20201104_at.least.one.is.AShM.beta0.pic/normal_down_tumor_down.intersect.txt",head=T,sep = "\t")
id2=read.table("../20201104_at.least.one.is.AShM.beta0.pic/normal_up_tumor_up.intersect.txt",head=T,sep = "\t")
id=c(id1$unitID,id2$unitID)

i=7
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
rt=file[file$unitID %in% id,]
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 8:10) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% id,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
file$mean.normal.beta0=rowMeans(file[,seq(2,16,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,16,4)],na.rm = T)
std <- function(x) sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))###算标准误差
file$std.normal.beta0=apply(file[,seq(2,16,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,16,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
original.normaldata=read.table("../20201104_at.least.one.is.AShM.beta0.pic/normaldata.txt",head=T,sep="\t")
orders=data.frame(unitID=original.normaldata$unitID,original.normaldata[,5:7])
file=merge(file,orders,by="unitID")

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)

p3=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""

p4=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")

layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p3,p4),layout=layout)
